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Abstract 

We investigate in detail the initial susceptibility, magnetization curves, and microstructure of 
ferrofluids in various concentration and particle dipole moment ranges by means of molecular 
dynamics simulations. We use the Ewald summation for the long-range dipolar interactions, take 
explicitly into account the translational and rotational degrees of freedom, coupled to a Langevin 
thermostat. When the dipolar interaction energy is comparable with the thermal energy, the 
simulation results on the magnetization properties agree with the theoretical predictions very well. 
For stronger dipolar couplings, however, we flnd systematic deviations from the theoretical curves. 
We analyze in detail the observed microstructure of the fluids under different conditions. The 
formation of clusters is found to enhance the magnetization at weak fields and thus leads to a 
larger initial susceptibility. The influence of the particle aggregation is isolated by studying ferro- 
solids, which consist of magnetic dipoles frozen in at random locations but which are free to rotate. 
Due to the artiflcial suppression of clusters in ferro-solids the observed susceptibility is considerably 
lowered when compared to ferrofluids. 

PACS numbers: PACS: 61.20.Ja Computer simulation of liquid structure; 75. 05. Mm Magnetic fluids; 
82.70.Dd Colloids; 82.20.Wt Computational modeling;simulation 
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I. INTRODUCTION 



Ferrofluids (dipolar magnetic fluids) are colloidal suspensions of ferromagnetic particles 
of about lOnm diameter dispersed in a carrier liquid that are usually stabilized against 
agglomeration by coating particles with long-chain molecules (sterically) or decorating them 
with charged groups (electrostatically). The small size of the particles favors magnetic mono- 
domains with a magnetic moment proportional to the volume of the magnetic grains. As 
a result, the particles interact with each other by the long-range anisotropic dipole-dipole 
potential as well as the short-range symmetric potentials, such as the steric repulsion, the 
electrostatic repulsion and the van der Waals attraction. The study of both the magnetiza- 
tion properties and structure of ferrofluids are of both fundamental and application interests. 

A ferrofluid system of sufficiently low concentration behaves like an ideal paramagnetic 
gas. The interactions between the particles can be neglected, and the physical properties 
of the system are well described by the one-particle model p|. In this case the equilibrium 
magnetization is written in terms of the Langevin function C{a) = coth(a) — 1/a 

N TJX 

Ml = MsatC{a), M^at = 77 — . (1) 

V fiQ 

Here, Mgat is the saturation magnetization of the fluid, a = mH/kT is the Langevin param- 
eter, and m is the magnetic moment, N/V denotes the number density of the particles, H is 
the magnetic fleld, T the temperature, k the Boltzmann constant, and fiQ = An x lO^^iJ/m. 
Eq.(|I]) leads to Curie's law for the initial susceptibility xl which depends linearly on the 
particle concentration 

_N rri^ 

~ vsil^' ^ ^ 

However, experiments on concentrated ferrofluids have revealed an essential deviation from 
the Langevin formula [§] . The initial susceptibility increases with the particle concentration 
much faster than that predicted in Eq.(||). This deviation is clearly due to the particle 
interactions. 

The influence of the particle interactions on the physical properties of ferrofluids is in- 
vestigated in terms of two dimensionless parameters. One is the volume fraction of the 
particles 
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with a the particle diameter. The other is the dipolar coupling constant 



2 

m 



which relates the dipole-dipole interaction energy of two contacting particles to the thermal 
energy kT. The two parameters can be combined into the Langevin initial susceptibility xl 
in Eq.(l) as Xl = 80 A. 

A number of theoretical models allow the evaluation of the equilibrium magnetization 
of ferrofiuids 0, ^, ||, ^, ||, |10|, |ll|. Among them, the mean-spherical model P, 0, the 



thermodynamic perturbation model and the modified variant of the effective field model 
have been shown to give good results for the magnetic properties of ferrofiuids with low 
or moderate concentration of magnetic particles (up to 10 — 12%). Taking xl as a universal 
parameter in weak magnetic fields, all the above approaches lead to the same expression for 
the initial susceptibility 

X = Xl(1 + Xl/3). (5) 

For ferrofiuids of higher concentration, larger particle size, or at low temperature, Eq.(^ is 
insufficient [0]. Higher order corrections in the dependence of x or M on (p, A, or in general 



on Xl, turn out to be crucial. 

An expression for x up to cubic accuracy in xl can be obtained by the mean-spherical 
model iP, the Born-Mayer expansion method [11 1 and the statistical model based on the 
pair correlation function [|l^, which reads as[|lO[ 

X = XL(l + XL/3 + xi/144). (6) 

The cubic term makes a noticeable difference of 5% over Eq.(|^) for xl > 4.14. The sig- 
nificance of the third order contribution has been verified in comparison with experimental 
measurements on the temperature dependence of x [0 • However, in what we classify here 
as the strong coupling limit, A > 2, the enhanced particle aggregation is not sufficiently 
accounted for by the analytical studies, thus giving rise to considerable deviations from the 
predictions of Eqs.(H,|D. Here the magnetic properties turn out to depend on $ and A sep- 
arately rather than xl alone, a point which we will try to elucidate with our simulation 
study. 

Previous simulation works on dipolar (ferro) fluids were mainly undertaken to explore 
the phase diagram. Large ranges of particle concentration and dipolar coupling strength 



have been examined. The results are found to depend substantially on the employed short 
range interactions, such as hard sphere |]12|, 0, |T^, soft sphere |T^, |1^, |T9[ and 
Lennard- Jones (LJ) potentials As an example, the question whether a minimum 

amount of dispersive energy, i. e., the attractive van der Waals energy, is required to observe 
the liquid- vapor coexistence is still under discussion |T^, |TB|, [TP|, Recent 
studies on the magnetic properties include the Monte Carlo simulations on the dispersions 
of interacting super-paramagnetic particles in a solid matrix and on ferrofluids with 
a finite spherical boundary [^. In order to compare with the macroscopic theories on 
ferrofluids, it is necessary to simulate the ferrofiuid systems with truly periodic boundary 
conditions. 

In this paper, the Langevin dynamics simulation method is used to study the equilib- 
rium properties of ferrofluids. The long-range dipolar interactions are dealt with the Ewald 
summation. The initial susceptibility and the magnetization curves are calculated as a func- 
tion of and A. The simulation results are directly compared with the theoretical models. 
The influence of particle aggregation is investigated systematically by performing a clus- 
ter analysis of the microstructure of the systems. The formation of aggregates is found to 
enhance the magnetization of ferrofluids at weak magnetic fields. This is confirmed by com- 
parison with the calculated initial susceptibility for a ferro-solid, i. e., a ferrofiuid with the 
same properties but with the particles' translational degree of freedom frozen in at random 
positions. 

The paper is organized as follows. We describe the simulation method in Sec. II. The 
simulation results on the initial susceptibility and the magnetization curves are given in 
Sec. Ill (A) and (B), respectively. The cluster analysis of the microstructure is presented in 
Sec. Ill (C). The initial susceptibility of the ferro-solid systems is calculated in Sec.III(D). 
We end with our conclusions in Sec. IV. 



II. SIMULATION METHOD 

The investigated ferrofiuid systems consist of spherical particles of diameter a dis- 
tributed in a cubic simulation box of side length L. Each particle has a permanent point 
dipole moment nii at its center. Using periodic boundary conditions in all spatial directions. 
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the dipole-dipole interaction potential between particle i and j is given by 

jrdip ^ 1 / ■ _ 3[mi ■ (r^j + nL)][mj ■ (r^^ + nL)] ) 

where rjj = — is the displacement vector of the two particles. The sum extends over all 
simple cubic lattice points, n = {n;j.,ny,nz) with n^, Uy, integers. 

In this work, we use the Ewald summation for dipolar systems to evaluate Eq.(|^) effec- 



tively, which gives [PGj ^ ^ 



= f/M + + U^f^^ + Ul^^^ , (8) 

where the real-space t/,^\the k-space (reciprocal-space) U^'^\ the self u'l^'^^^\ and the surface 
jjipi~f) contributions are respectively given by: 

Ui;^ = I] { (m^ • m,)5(| r,, + n I) - [m, ■ (r,, + n)] [m, ■ (r,, + n)]C(| r,, + n |)(l 

^^^^ " 47r/!oL3 5Z ^exp[-(7rA;/KL)2](mi ■k)(mj- ■k)exp(27rik-rij/L), 



(10) 



keZ3,k7^0 



f/,^r^^ = ^7 , ■ m,. (12) 



with 

B{r) = [erfc(Kr) + (2/^r/v^)exp(-/^V)]/r^ (13) 
^(r) = [3erfc(Kr) + (2Kr/v^)(3 + 2KV2)exp(-KV)]/r^ (14) 

Here erfc(x) := 27r~^/^ JJ^ exp(— t^)dt denotes the complementary error function. The in- 
verse length K is the splitting parameter of the Ewald summation. Eq. (|l^) assumes that the 
set of the periodic replications of the simulation box tends in a spherical way towards an 
infinite cluster and that the medium outside this sphere is a uniform medium with magnetic 
permeability fiBC- In this work, we use the metallic boundary condition fiBc = oo. The sur- 
face term vanishes in this case and demagnetization effects do not occur. Thus the applied 
external magnetic field H coincides exactly with the internal field. The influence of different 
boundary conditions on the simulation results will be the topic of a subsequent publication 
p9| . The related formulas for the dipolar forces and torques can be found in a previous 
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paper |2^. The theoretical estimates of the cutoff errors in the Ewald summation derived 
there are used to determine the optimal values for the Ewald parameters which enable us 
to minimize the overall computational time at a predefined accuracy. 

The short range interaction potential between the particles is chosen to mimic the steric 
mechanism for stabilizing the colloidal suspension. As mentioned in Sec. I, it is still not 
settled which kind of potential is most suitable to model real ferrofluids. It might actually 
depend on whether the solution is electrically or sterically stabilized. In this work we adopt 
the Lennard- Jones potential 

12 / _ \ 6 



LJ 



4e 



(15) 

where C{Rc) = {(t/RcY^ - {(y/Rcf with a cutoff radius of R^ = 2^/6^. In this way the 
particles have a purely repulsive interaction force which smoothly decays to zero at Re- 
This cutoff range is smaller than the one usually used in the soft sphere potential which is 
commonly adopted as Uss{fij 



Ae{a/rijf'^ 1T|, [T^, |T8|, |T9[, and thus closer to the hard 



sphere potential as employed in most other theoretical calculations |T0|, ITl|. In this context 



the analytical calculations of Ivanov fTO] reveal that at least for the spatially homogeneous 
systems the influence of different short range repulsive potentials seems to be marginal. 
The translational and rotational Langevin equations of motion of particle i are given by 



27, 30 



I,; ■ ij^i 



T 

i 1 



R 
i 1 



(16) 
(17) 



where M.i and I, are the mass and inertia tensor of the particle, Tt and Tr are the transla- 
tional and rotational friction constants, respectively. The Gaussian random force and torque 
have the properties that the first moments vanish 



< il{t) > = 0, 



0, 



(18) 
(19) 



while the second moments satisfy 



(20) 
(21) 



where a and (3 denote the x, y, z components in the Cartesian coordinates, respectively. 
Introducing the dipolar and short range terms into Eqs.(|l^) and (p!7[), the dimensionless 
equations of motion can now be written as: 

V* = + ^f/*) - TtV* + (22) 

n-^: = Yl + < X H* - V\u^^ + if\ (23) 

where the variables are given in dimensionless form reduced by the following units: length 
r* = r/cr, dipole moment m*^ = m^/47r/ioe(T^, moment of inertia /* = //(A^cr^), time 
r = t(e/A<(T2)i/2, the friction constants = Yria'^ j Mefl^ and Y\ = TR/{M(x^ey/^, 
magnetic field H* = H^AnfiQcr^ /eY^'^ as well as temperature T* = kT/e. The values of the 
dimensionless friction constants do not affect the equilibrium properties. Here we adopt 
r^- = 10.0 and T*^ = 3.0, because in our simulations these values had been found to speed 
up the equilibration time. These values are also close to that used in Ref. [^|. A value of 



/* = 0.4 is used for the dimensionless moment of inertia, corresponding to that of a rigid 
sphere with diameter a. 

The simulations were performed at constant temperature T* = 1. The orientational coor- 
dinates of the particles were expressed in terms of quaternion parameters and the equations 



of motions were integrated with a leap-frog algorithm [^. A reduced time step At* = 0.002 
was employed in all simulations. The runs were started from initial configurations with 
random particle positions and dipole moment orientations. For each case, the system was 
at first equilibrated for a period of 50 000 time steps. The magnetization and structural 
properties were then calculated from the data for another period of at least 200 000 time 
steps. Error bars for the simulation results were determined by dividing the simulation runs 
into blocks of 10 000 time steps and calculating an estimate for the standard deviation of 



the mean 31 . 



III. RESULTS AND DISCUSSION 



Results were obtained from simulations on systems with N = 1000 particles. The number 
density of the particles was mostly taken in the range of p* = N/V* < 0.45, corresponding 
to volume fractions (p = p*7r/6 < 0.236. This covers the typical particle concentration ranges 
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which are obtainable in real ferrofluids. A few more runs were also carried out at higher 
concentrations in case of smaller A in order to reach a sufficient large xl- The values of 
the dipolar coupling constant A were mainly chosen to be 1,2,3 and 4. To get an estimate 
of realistic A-values, we consider a ferrofluid sample consisting of monodispersed magnetite 
particles. The bulk magnetization of this material is = 4.46 x 10^ A/m. With a magnetic 
core diameter of the particles of cr = 10 nm, the value of A amounts to about 1.3 at room 
temperature T = 3007^^. For larger particles with a = 13 nm, A rises up to 2.9 under the 
same conditions. Note that A may further be enlarged by decreasing the temperature or 
using cobalt as working material (M^ = 14.01 x 10^ A/m), for which the previous A-values 
increase by a factor of around 3.14. We therefore also carried out some simulations for A = 8. 

In all simulations, we fixed the root mean square (rms) absolute errors in the dipolar 
forces to AF'^^^ < 10~^m^/47r/ioO"^, which corresponds to 10~^/6 of the attractive force 
between two contacting particles with dipole moments in parallel alignment. The optimal 
values of the Ewald parameters had been determined separately for each system . As an 
example, for a sample with p* = 0.1 we get Vc = 0.425L, k = 6.78/ L and Kc = Q x 2%/ L. 
While p* = 0.3 yields Vc = 0.46L, k = 6.84/ L and Kc = 7 x 2n/L to be the optimum values, 
respectively. 



A. Initial susceptibility 

The infiuence of the particle interactions on the magnetic properties is most evident at 
weak fields, so we start with the investigation of the initial susceptibility x- The value of 
X is determined from the linear relationship of M = aX H ^ 0. Simulations are thus 
performed to get the M{H) curves at weak fields. The equilibrium magnetization M is 
evaluated by the prescription 



M = -/ (— -^m,)rft. 



(24) 



where to is the time needed for the system to establish the equilibrium state. 

To study the weak field magnetization, we increase the Langevin parameter a from to 
0.5. As expected, the calculation of M at a = does not show the existence of a spontaneous 



magnetization, because is still not high enough. Other works [19| have shown that the 



transition to a ferromagnetic state occurs at ~ 0.44 for A = 4 . For higher values of A, the 
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critical concentration is around ~ 0.35 for A = 6.25 |T3|, [T^ and ~ 0.31 for A = 9 [|T7 . 
These conditions are clearly not encountered in our simulations. 

The dimensionless magnetization M*(= M/ Ane/ ^^a^) of the system with A = 3 is 
shown in Fig.l as a function of a and 0. The results for other values of A display a similar 
behavior. It can be seen that M* at first increases linearly with a, then at higher field 
values sub-linearly. We introduce aun as a definition for the maximum value of a where M* 
still can be considered to obey the linear relationship. Evidently, at a fixed value of A this 
a^jn- value decreases with growing 0. The data in Fig.l show that when is increased from 
0.052 to 0.209, aiin decreases from 0.45 to about 0.25. On the other hand, at a fixed particle 
concentration, rising A also implies a reduction of aun- This demonstrates that the range 
over which M is proportional to H (or equivalently a) shrinks with the increase of xl- For 
a quantitative evaluation of the initial susceptibility we need to determine the linear region 
for each case separately. This is accomplished with the help of a linear regression fit with 



error weighting The basic idea of that method is to minimize the so-called ^-square 
merit function 

2 ^ M.-M(0)-xg. .2 
1=1 

so as to get x ^"^^ its error bar. Here AMj is the error bar related to the simulation data 
of Mj. In Fig.l the data were taken with a step size of Aa = 0.05 giving 11 data points 
on each curve. Since no spontaneous magnetization occurs, the initial data point is always 
at M(0) = 0.0. For each curve, a series of linear fittings was generated by increasing in 
Eq.(p5D from 5 to 11. The fit which gives the minimum value of per degree of freedom 
is taken to be the optimum linear fitting to the data. The associated A^^^-value is referred 
as Nd,opt- The slope x aiid its error bar determine the initial susceptibility. The value 
of aiin is roughly given by the product of Aa with N^^opt — 1- The results of this data fit 
including the error estimates are compiled in Table. I. The degree of reliability of the results 
is reflected by the high goodness-probability Q [^. The resulting straight lines M = xH, 
or in dimensionless form M* = xa/47r a/X/T* , are also displayed in Fig.l to show how they 
fit to the simulation data. 

Fig. 2 (a) shows the calculated initial susceptibility x as a function of A and 0, along with 
the analytical results from Eqs.(|]) and (Q). Since the Langevin susceptibility xl has been 
taken as the universal parameter in these theoretical predictions P, 0, 110], we also map the 



data of X oii Xl in Fig. 2(b) to show the scahng behavior more directly. In addition to the 
results in Fig. 2(a), some more data was produced at higher concentrations (up to = 0.419) 
for the case of A = 1 in order to get a larger range of xl- We did not combine such high 
concentrations with larger values of A, since the close proximity to the magnetic phase 
separation would make the discussion more complex, rendering the results more sensitive to 
the specific short-range interaction. 

For the weak interaction range A < 2 we find perfect agreement with Eqs.(|D and (^ 
for all investigated volume fractions. Actually, because for these coupling parameters we 
can reach only up to xl < 4, both analytical predictions yield equally good results for our 
data and are indistinguishable within the error bars. For A = 3 we see small but significant 
deviations at low values of 0, meaning that already at low values of xl neither Eqs.(^ nor 
(D are accurate. However, at high volume fractions (xl > 3), the third order equation (H) 
is compatible within error bars. For the numerical values of the data see Table 1, where also 
the values of the analytical results for different orders are given for comparison. Even more 
drastic deviations from the theoretical predictions are found for A = 4 at low concentrations, 
in the range Xl ^ 6. We observe again the trend, that at higher volume fractions the third 
order equation (|^) becomes compatible with our data. We can conclude that the parameter 
A, rather than xl-, becomes the dominating control parameter for A > 2. This effect can be 
attributed directly to the microstructure formation in the system. When A < 2, what we 
consider to be the weak interaction limit, the dipolar potential is comparable to the thermal 
fluctuation. Then the interparticle correlations are weak, and the present theoretical models 
work well. However, going to stronger interactions with A > 2, leads to a considerable 
amount of particle aggregation. The resulting clusters amplify the sensitivity of the system 
to weak external magnetic fields and thus lead to an enhanced magnetic susceptibility. We 
will come back to this point in Sec. pi C. 



The influence of the volume fraction seems to be weak in the investigated parameter 
regime. Noteworthy is that the third order equation seems to describe the available data 
well at larger volume fraction and high coupling parameter A, even better than it does at 
intermediate values. On the basis of our available data we cannot decide if this finding 
is a mere coincidence. To our knowledge there is no theoretical reason to believe that 
the third order equation should become reliable again at these high values of xli where 
higher order corrections of the form 0'^A" with m, n arbitrary, are expected to become 
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important. Corrections linear in 0, but up to ninth order in A have been assembled for 
the Born-Mayer expansion method in Ref. |Tl|. These additional terms are all positive 
and therefore necessarily increase the susceptibility. We have found that they improve the 
theoretical prediction of x only at very small concentrations {(f) < 0.05) for A > 2, but they 
dramatically overestimate the simulation results at larger values of 0. 

Note that the simulations presented in Ref give rise to a susceptibility x that decreases 



with growing A if xl is fixed, which is opposite to our findings in Fig.2(b). We suspect that 
this trend may be induced by finite size effects, since a finite spherical simulation volume 
with only a few hundred particles had been used in their simulations (without a periodic 
replication to properly account for the long range dipolar interaction as is done in the present 
work). We shall clarify this point in a subsequent publication ||29||. 



B. Magnetization curves 

In this section, we extend the calculation of M to a larger range of the applied external 
field strength a in order to get the full magnetization curves. For moderately concentrated 
ferrofluids, the 'modified mean-field model' |^ describes the magnetization curve in terms 
of the Langevin magnetization Ml [Eq.(|l])] as follows 

M{H) = Ml{H,), (26) 

where the effective field acting on an individual ferro-particle is given by 

H, = H + Ml{H)/?>, (27) 

rather than = H + M/3 as done in the Weiss model. The resulting expression for the 
initial susceptibility is Eq.(|). A higher order description is obtained by replacing He in 
Eq.(gZl) by 

H^ = H+'Ji^ + ^MUH)'-Mi^, (28) 



which in the weak field limit reduces to Eq. (y) for x 1131 • '^^^ effective Langevin parameters 



corresponding to Eqs.(p^) and (|28|) can be written as = a + XL^ict) and = a + 



XL>C(a;) + ^£(a)^^^, respectively. This reflects the x^-dependence of these two theoretical 
predictions. The Born-Mayer expansion described x^^ ^ different way, where its dependence 
on and A are given separately [pJ]] . 
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We investigated the magnetization curves for two different cases. In the first case xl 
is fixed to a small value of 1.256 where the predictions using Eqs.(p7D and (|28| ) practically 
coincide. Then A is taken to be 1,2,4 and 8, respectively. To keep xl constant, the volume 
fraction (p had to be adapted from 0.157 down to 0.0196. In the second case, we take a larger 
value of Xl = 5.026 which makes the difference between the two theoretical predictions 
distinguishable by a few percentage. Only A = 3 and 4 are discussed. The corresponding 
volume fractions are (p = 0.209 and 0.157, respectively. In the latter case, we did not include 
the results for A < 2 due to the limitation of the particle concentrations and for A = 8 due 
to the extremely slow dynamics which prohibited us from reaching equilibrium distributions 
within our computational limits. 

The calculated magnetization curves are shown in Fig. 3 together with the theoretical 



predictions using Eqs.( P6H28| ). For a more explicit demonstration of the particle interaction 
the inlays display the difference between the simulation data and the Langevin model (M — 
Mij/Msat- The characteristic maximum of this difference is found to be at intermediate field 
strength of a ~ 1, where it may reach more than 35%. The results in Fig. 3 (a) reveal again 
the essential role of the dipolar coupling constant A. When A equals 1 and 2, the simulation 
and theoretical curves agree with each other very well. But for larger A the calculated 
magnetization M rises with a much faster than the theoretical predictions, especially in the 
weak field regime. This reflects the same result already observed in Fig.2, that x increases 
with A at low and moderate concentrations if xl is held fixed. Due to this effect, the peak 
positions of the increment {M — Ml) / Mgat shift to smaller values of a with the increase of A. 
In Fig.3(b), the volume fraction of the particles are relatively high. The maximum absolute 
difference between the two analytical predictions occurs at a = 0.62 where the result given 
by using Eq . (|28D is about 4.3% larger than that by Eq . (|27|) . Our simulation results are found 
to be consistent with both of them within the range of error bar. A quantitative analysis 
of the data indicates that Eq. (p8D gives a better description in the weak field regime. This 
consequently leads to a better agreement between the simulation results on x with the 
prediction of Eq.(^. However, both data sets for A = 3 and A = 4 show deviations from 
the theoretical curves at larger values of a which cannot be explained by statistical errors. 
Also here we infer that higher order corrections are needed to more accurately describe the 
magnetization curves. 
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C. Microstructure analysis 



The microstructure and its relation to the magnetic properties of ferrofluids are investi- 
gated by performing a cluster analysis. In conventional cluster analysis, agglomerates are 
usually defined on the base of the spatial proximity between the particles or by means of 
an energy criterion. For ferrofluids the latter one is more favorable, because the anisotropic 
dipolar interaction implies that two neighboring particles can form a stable bonding only if 
their dipole moments are roughly aligned in head-to-tail orientation. 

The energetically based cluster analysis of ferrofluids has been done in several different 
ways n, TU, 33 1. In this work, we adopt the definition of Refs. |T2|, H and W 



where two particles are considered to be bound if their dipolar potential energy is less than 
a predetermined value Ubond- The influence of using different threshold values for Ubond has 



been investigated in a recent work [Q. The values of — 1.4AA;T |p!2| , and — 1.5AA;T |19 
were shown both to give good results. In our calculations we have varied Ubond from —lAXkT 
to —l.QXkT. Only slight quantitative changes of the results have been found without a 
quahtative difference. The results shown here have been evaluated for Ubond = — l.SAfcT, 
i. e., 75% of the contact energy of two perfectly co-aligned dipolar particles. The above 
cluster definition implies that aggregates are mainly linear chains and rings. Thicker particle 
agglomerates or branched chains, however, occur very rarely in our simulations. 

The size of a cluster is defined as the number of particles belonging to it. We first set up 
the connectivity matrix by means of our cluster definition, and then use it to evaluate the 
sizes and numbers of clusters in a given configuration. The average cluster size is defined by 

ris > ■ (29) 

s s 

where Ug is the number of clusters having size s, and the triangular brackets denote the time 
average, or equivalently, the average over the configuration space. 

Fig. 4 plots Savg as a function of A and for the case of zero field {a = 0). Clearly, 
increasing the dipolar coupling A leads to a larger cluster size Savg- However, in the weak 
coupling limit A = 1 or 2, Savg exceeds 1 only slightly, meaning that there are very few 
dimers present. Inter-particle correlations seem to be weak and we already noted that the 
simulation results agree well with the theoretical models in these cases. 

As A is increased, the larger values of Savg indicate the formation of more clusters. To 
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elucidate this, we show snapshots of the three-dimensional configurations in Fig. 5 for the case 
of = 0.052 (p* = 0.1) at different A values. The increase of the number and size of clusters 
with A can be seen clearly in these figures. For A = 3 we find Savg = 1-123 at = 0.052 
. About 10% of the particles are organized in dimers and less than 1% in trimers. At the 
same density, Savg increases to 1.481 for A = 4. Correspondingly, 22% of the particles are in 
dimers, 7% in trimers and 2% in quadrumers. The density dependence of these numbers is 
shown in Fig. 6 where the average percentage of particles in ra-mers are plotted as a function 
of for A = 3 and 4. It can be seen that more than 60% of the particles are still in monomers 
and the clusters with s > 5 are rarely found in these cases. However, when A equals to 8, all 
the particles get involved in the formation of long chain-like structures. The average chain 
length is larger than 13o" at = 0.052. Due to thermal fiuctuations, these chains are quite 



fiexible. They are continuously breaking and recombining with others |]T2|, |14|. Even closed 
ring structures can be observed occasionally. In this strong coupling case, the chains have a 



close analogy to living polymers M, 34 



When the external field is switched on, the tendency of the dipole moments to co-align 
with the field direction enhances the aggregation probability. Accordingly the average cluster 
size is expected to increase with the field strength. Fig. 7 shows the field-dependence of Savg 
for the systems studied in Fig.3(a). The increment of Savg with a is only very weak for 
A < 2, but it is more evident at A = 4. As exhibited in Fig. 8, the structure character of 
the system varies from randomly distributed short clusters to string-like alignments along 
the magnetic fiux with the increase of a. Consequently the cluster analysis in Fig. 9 refiects 
that more and more particles merge into larger aggregates. Cluster as long as s = 15 could 
be observed at high a. But these string-like structures in Fig. 8(d) are much more unstable 
than the chains formed at A = 8 [see Fig.5(d)]. They posses a smaller contact energy so 
that thermal fiuctuations can easily break long clusters into smaller segments. This effect 
is field- independent. In fact, most of the strings observed in Fig. 8(d) consist of a few 
continuously breaking and recombining segments. As a result, the average cluster size Savg 
does not present a dramatic change as expected from the snapshot of the configurations. 
When A = 8, however, the long-chains had already been formed even at zero field. The 
increase of the field strength just arrange these chains to have a better alignment along the 
field direction. But the average chain length does not change too much. 

The results in Fig. 7 reveal that the magnetic field dependence of Savg remains weak at 
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small a. It is therefore legal to invoke the zero- field cluster analysis in order to interpret 
the outcomes related to the initial susceptibility (Fig. 2). Let us focus to the cases of A = 3 
and 4 at low concentrations where the deviation between theory and simulation are most 
pronounced. In these cases the small clusters are relatively stable and spatially well separated 
from each other. For a rough approximation, we treat the clusters as rigid particles thus 
considering our system as a poly-dispersed sample. We define the effective dipole moment 
of each cluster according to 



s 

i=\ 



= 5^m,. (30) 

The effective Langevin susceptibility xl^ of such a poly-dispersed system can be computed 



as 



Since most of the clusters are short linear n-mers, the value of Xl"^ larger than xl - 
the appropriate value for the mono-dispersed system. Replacing xl by x'l'^ Eq.(^ we 
get new modified theoretical curves. In Fig. 10 it can be seen that this simple operating 
prescription considerably improves the agreement between simulation and theory, at least 
at low concentrations. This proves that the formation of small clusters is responsible for the 
enhancement of the initial susceptibility. 

At particle concentration above ~ 0.1, the modified theoretical curves start to overesti- 
mate the simulated x-values. In these cases, the increased angular correlations between the 
particles are not just limited to the two next nearest neighbors and thus make the above 



definition of clusters ambiguous [0. Although we keep on using the energy criterion to 
obtain the cluster informations for Fig.4 and 6, these values are more likely refiecting the 
local orientational order of the dipole moments. The clusters we find are not stable entities 
involving well-defined particle numbers as in the low density regime. This can be partially 
seen from the Savg{4>) behavior in Fig.4. When A < 2, Savg slightly increases with 0. But 
it saturates at A = 4 and even decreases at A = 8 for > 0.1. Next we examine the radial 
distribution function g{r) of the system (Fig. 11), which gives the probability of finding a 
pair of particles a distance r apart, relative to the probability expected for a completely 
random distribution at the same density [^. We find that the height of the first and second 
peaks, respectively, associated with dimer and trimers, decrease with the increase of 0. This 
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implies that the number of clusters does not increase as rapidly as the density |T9[ . In other 
words, upon rising (p the spatial distribution of the particle positions re-approaches to that 
of a homogeneous system. Therefore we expect our suggested modification with Eq.(pTD 
only to be valid at low volume fractions below cj) = 0.1. 

The formation of long chain-like structures in the strong coupling regime makes the 
systematic investigation of the magnetization properties at A = 8 difficult. At = 0.0196, 
there are very few overlaps between the elongated chains. They have enough free volume to 
re-arrange in response to the external field. The magnetization of the system increases very 
fast in the weak field regime, as shown in Fig.3(a). The occasional formation of a few closed 
ring structures is not found to affect noticeably the global magnetization behavior. On the 
other hand, when is enlarged, the chains become strongly entangled, like in a semi-flexible 
polymer solution. This slows down the dynamics of the system dramatically. An extremely 
long simulation time is required to reach the equilibrium state in the weak fleld regime. For 
example, at a density of = 0.0785, the equilibrium value of M at a = 0.5 is obtained only 
after 1.5 million simulation time steps. Simultaneously, since the number of particles is fixed 
in our simulation, the increase of makes the chain length become comparable to the size of 
the simulation box. This unavoidably induces finite size effects. For these reasons, we were 
unable to provide systematic investigations at A = 8 at higher concentrations (or equivalently 
larger Xl), but only for xl = 1-256 shown in Fig.3(a). These numerical difficulties could 
well be related to experimental observations which find a dramatic increase of viscosities at 
similar dipolar couphngs p5| . 



D. Comparison with ferro-solids 

To further investigate the infiuence of the cluster formation process we now turn to ferro- 
solid systems. By this we denote a suspension of dipolar particles which have their spatial 
positions frozen in, but are permitted to rotate freely. In this way the aggregation related 
effects in ferrofluids can easily be separated. 

We prepared the initial conflgurations of the ferro-solid system by placing the particles 
randomly in the simulation box with random dipole moment orientations. For each given 
A and 0, several different initial conflgurations were generated and simulated. They were 
found to give very similar results. Computer time limitations prohibited us to average over 
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a large sequence of quenched configurations, which would be desirable to obtain a good 
ensemble average for a ferro-glass system. However, we have investigated this case only to 
study the qualitative influence of freezing the positional degrees of freedom. Accordingly, 
only the rotational equations of motion [Eq.(|T^] of the particles need to be integrated in the 
simulations. Fig. 12 displays the simulation results for the initial susceptibility at A = 3 and 
4, along with the theoretical predictions for the fluid system. In ferro-solids x ^i-lso exhibits 
a nonlinear increase with the volume fraction cf). Comparing our data with the results of the 
ferrofiuid system in Fig. 2, we see that the suppression of clusters leads to a slower increase 
of the magnetization at weak fields, and hence to a lower x- K is interesting to note that the 
discrepancy between XsoUd and x fluid goes up when A or are increased. This shows that 
the positional correlations are also very important at high density, and couple strongly to 
the dipolar rotational degrees of freedom, even if they do not evidently form stable clusters. 
We consider the fact that the measured values of XsoUd coincide with the theoretical curves 
for the fluid system, Eqs.(^ and @, almost exactly at low concentration for an accidental 
correspondence. At higher volume fractions the ferro-solid displays a susceptibility even 
below the theoretical prediction. This effect is more evident for larger A at higher 0. Note, 



that in Ref. ^ a similar solidified ferro-colloid was investigated. There XsoUd was found to 
increase with xl-, but the values are larger than that for the associated ferrofiuid. Again we 
suspect that this due to finite size effects of their sample. 

We point out that the particles in our ferro-solid systems are supposed to be super- 
paramagnetic. This means that they have a negligible crystallographic magnetic anisotropy 
energy. There is no energy barrier which prevents the magnetic dipoles from a free rotation. 



That makes the present systems different from that one studied in Ref.|2J where the initial 
susceptibility is found to decrease with the increase of the particle concentrations. 



IV. CONCLUSION 



We have investigated in detail the initial susceptibility, magnetization curve and mi- 
crostructure of ferrofluids in various particle concentration and dipolar coupling strength 
ranges by means of molecular dynamics simulations. The ferrofiuid was modeled as a soft 
sphere system by means of a purely repulsive Lennard- Jones potential, and the magnetic 
cores interacting with a dipolar coupling strength A from 1 to 8. The investigated volume 
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fractions were in the range of < 0.24 mostly. Our simulation method used the Ewald 
summation with metallic boundary conditions to deal with the long-range dipolar interac- 
tions, and took explicitly into account the translational and rotational degrees of freedom 
of the dipolar particles. The temperature was kept constant by means of a Langevin ther- 
mostat for all degrees of freedom. Our simulation results for the susceptibility show good 
agreement with the theoretical predictions in the weak coupling limit A < 2. For higher 
dipolar couplings, however, we find systematic deviations from the theoretical predictions at 
low and intermediate concentrations. The detailed cluster analysis demonstrates that this 
is due to the particle aggregation at higher A. The formation of clusters tend to increase the 
magnetization at weak field regime and consequently induce a larger initial susceptibility. At 
high densities, the spatial distribution of the particles starts to homogenize again, and the 
significance of clusters goes down. In fact we speculate that at higher concentrations there 
are almost no stable clusters any more, because they can disintegrate and reform with other 
particles rapidly. In this regime we find that the third order equation in xl is compatible 
with our data. 

For the magnetization curves we find that both, the second and third order theoretical 
curves, describe our data well for couplings A < 2. For higher values of A we find again 
systematic deviations. The system shows at coupling strength of A = 8 an almost gel-like 
behavior due to entangled chain structures. We expect this effect to become experimentally 
relevant in Cobalt based ferrofluids, where A is large. 

The influence of the aggregation phenomena was isolated by studying ferro-solids where 
the particles are randomly placed inside the simulation volume, and only rotations are 
allowed to occur. The observed initial susceptibilities are lower compared to the ferrofiuid 
systems due to the suppression of cluster formation. The deviations to the fluid system grow 
with increasing dipolar coupling strength and increasing volume fraction. 
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TABLE I: Linear regression fitting results for the M(H) [M(a)] curves in Fig.l. The theoretical 
predictions of x o^i different orders are also given for comparison. 
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0.963 


6 
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FIG. 1: Dimensionless magnetization M* as a function of the Langevin parameter a and the 
volume fraction cf) for a system with A = 3. The symbols are the simulation results. The dashed 
curves are the linear fitting data, the slope of which gives the initial susceptibility %. 
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FIG. 2: (a) Initial susceptibility x ^-s a function of the dipolar coupling constant A and volume 
fraction 0; (b) same data of x iii (3-)! but mapped on the Langevin susceptibility xl- The 
theoretical curves are given by the Langevin model (x = XL) dashed), Eq.(|5[) (dotted) and Eq.@ 
(dotted-dashed), respectively. 
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FIG. 3: Magnetization curves in case of (a) xl = 1-256 for A = 1 (circle), 2 (square), 4 (diamond) 
and 8 (triangle- up) ; (b) xl = 5.026 for A = 3 (circle) and 4 (diamond), respectively. The theoretical 
curves are given by the Langevin model [Eq.(|l|), dashed line], and by using Eq.(|^) (dotted line) 
and Eq.(p8|) (dotted-dashed line) in Eq. (|26|) , respectively. The inlays show the difference between 
the calculated magnetization and the Langevin model (M — ML)/Msat- 
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FIG. 4: Average size of the clusters formed at zero field as a function of A and 
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FIG. 5: Snapshot of the three-dimensional configurations formed at zero-field in the systems of 
4> = 0.052 with A = 1 (a), 3 (b), 4 (c) and 8 (d), respectively. The arrows with length of a indicate 
the orientation of the dipole moments. The spheres (diameter=cr/2) show the center positions of 
the particles. 
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FIG. 6: Average percentage of the particles in n-mers for the system with A = 3 (a) and 4 (b) at 
zero-field. 
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FIG. 7: Ficld-dcpcndence of the average cluster size of the systems studied in Fig.3(a): A = 1 
(circle), 2 (square), 4 (diamond) and 8 (triangle- up) . 
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FIG. 8: Snapshot of the three-dimensional configurations formed in the system with (j) = 0.0392 
and A = 4 [studied in Fig.3(a)] under the magnetic field of a = (a), 1.0 (b), 3.0 (c) and 9.0 (d), 
respectively. The arrows on the frame indicate the direction of the magnetic field. 
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FIG. 9: Field-dependence of the average percentage of particles in n-mers for the system with 
= 0.0392 and A = 4. 
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FIG. 10: Comparison of the simulation results on the initial susceptibility with the modified 
theoretical model. The curves are the theoretical prediction of Eq.(^) (dotted) and the modified 
theoretical results obtained by replacing xl with the effective Langevin susceptibility x^^^ in 
Eq. ^ (dotted-dashed) , respectively. 
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FIG. 12: Simulation results on the initial susceptibility of the ferro-solid systems with A = 3 and 4 
as a function of the volume fraction (j) and the dipolar coupling parameter A. The analytical results 
on the ferrofluid system from Eq.(^) (dotted) and Eq.(P) (dotted-dashed),as well as the simulation 
results from Fig. 2 are also included for comparison. 
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